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Abstract 

We present a numerical study of the flow of an as- 
sembly of frictionless, soft discs at zero tempera- 
ture, in the vicinity of and slightly above the jam- 
ming density. We find that some of the flow prop- 
erties, such as the fluctuations in the number of 
contacts or the shear modulus, display a critical 
like behaviour that is governed by the proximity to 
the jamming point. Dynamical correlations dur- 
ing a quasistatic deformation, however, are non 
critical and dominated by system size. At finite 
strain rates, these dynamical correlations acquire 
a finite, strain-rate dependent amplitude, that de- 
creases when approaching the jamming point from 
above. 
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1 Introduction 
ground 



and back- 



The properties of granular materials, in particu- 
lar sand, are a constant source of fascination for 
children and adults alike, and are intrinsically re- 
lated to the ability of such systems to exist in cither 
solid or fluid states, under very similar conditions. 
For the scientist this fascination may arise through 
the apparent contradiction between the rigidity of 
the individual grains and the fragility of the as- 
sembly as a whole. This means that, for exam- 
ple, small changes in the loading conditions (such 
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as changing the inclination angle of the support) 
can lead to large scale structural rearrangements 
( "avalanches" ) or even to the complete fluidization 
of the material. A few years ago, Liu and Nagel — 
have suggested a "phase-diagram" for this type of 
solid- liquid transition ("jamming"). At zero tem- 
perature the axes relate to the ways an unjamming 
transition can be triggered, either by increasing the 
external driving (e.g. shear stress) or by decreasing 
the density of the material. The present study will 
probe the vicinity of this "unjamming line" using 
quasistatic and finite strain rate simulations of a 
model granular system. 

If one applies a shear stress, which is small and 
below a certain threshold ( "yield-stress" ) , the ma- 
terial will respond as an elastic solid. Increasing the 
stress above the yield-stress, the particles will un- 
jam and start to flow. This flow behavior is called 
"plastic flow" as the material will not revert to its 
original shape when the stress is removed. 

In the following we will assume that the sys- 
tem can flow at arbitrary small strain-rates without 
showing flow-localization. That this is possible is 
by no means guaranteed, as in some instances co- 
existence of flowing and jammed states is observed 
£ We have not observed such a persistent strain 
localization in the simulations that are presented 
here. 

Plastic flow is observed in a large number of 
glassy materials, that are a priori very different 
from the athcrmal granular systems close to point J 
(see below) studied in this paper. However, all ma- 
terials display a rather universal behavior, that was 
illustrated already in early studies on the plastic- 
flow of metallic glasses .^4 These studies have given 
indications that in the flowing phase the main plas- 
tic activity is spatially localized to so called shear 
transformation zones.— £ These zones are non per- 
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sistcnt, localized in space, and presumably consist 
of a few atoms that undergo the irreversible re- 
arrangements responsible for the observed plastic 
flow. Recently, this plasticity has been further an- 
alyzed in simulations with a focus on the quasistatic 
dynamics at small strain rates, close to the flow ar- 
rest.— — With these studies it was possible to trace 
back the origin of plastic activity to the softening 
of a vibrational mode and the vanishing of the as- 
sociated frequency.— In real space, this softening 
is associated with the formation of distinct, local- 
ized zones where the plastic failure is nucleated.— 
In turn, this can trigger the failure of nearby zones, 
such that avalanches of plastic activity form that 
may "propagate" through the entire system. 

It has been argued that the macroscopic extent 
of these avalanches is a signature of the quasistatic 
dynamics, which gives the system enough time to 
propagate the failure throughout the system. Be- 
yond the quasistatic regime, i.e. farther away from 
the jammed state, the size of these events is ex- 
pected to be finite. Thus, one naturally finds an 
increasing length-scale that is connected with the 
flow arrest upon reducing the stress towards the 
threshold value.— ~— 

Without external drive, an (un)jamming transi- 
tion can occur for decreasing particle volume frac- 
tion below a critical value, <p c . This special point, 
which is only present in systems with purely repul- 
sive steric interactions, has been given the name 
"point J".i&±Z. At this point the average number 
of particle contacts jumps from a finite value zq to 
zero just below the transition. The value of zq is 
given by Maxwell's estimate for the rigidity transi- 
tion—^ and signals the fact that at point J each 
particle has just enough contacts for a rigid/solid 
state to exist. This marginally rigid state is called 
"isostatic" . Compressing the system above its iso- 
static state a number of non-trivial scaling prop- 
erties emerge. 16 ' 20 As the volume fraction is in- 
creased, additional contacts are generated accord- 
ing to Sz ~ The shear modululs scales as 
G <~ p/Sz and vanishes at the transition (unlike 
the bulk modulus).— This scaling is seen to be a 
consequence of the non-affine deformation response 
of the system,— with particles preferring to rotate 
around rather than to press into each other.— As- 
sociated with the breakdown of rigidity at point J 
is the length-scale, I* ~ fe -1 , 23 ' 24 which quantifies 
the size over which additional contacts stabilize the 



marginally rigid isostatic state. 

In this article we present results from quasistatic 
and small strain-rate flow simulations of a two- 
dimensional system in the vicinity of point J. To- 
gether with the linear elastic shear modulus, at 
point J also the yield-stress a y vanishes.— ~— Thus, 
point J is connected with a transition from plastic- 
flow behavior (<fi > <fi c , a v > 0) to normal fluid 
flow (<fi < cj) c , (jy = 0), with either Newtonian— 
or Bagnold rheolog y 27 ' 29 at small strain-rates. In 
consequence, both (un)jamming mechanisms as de- 
scribed above are present at the same time: the 
flow arrest, as experienced by lowering the stress 
towards threshold, is combined with the vanishing 
of the threshold itself. 

In this study we want to address two questions: 
In how far do the general plastic flow properties 
carry over to this situation of small or, indeed, van- 
ishing yield-stress? Is the vicinity to point J and its 
isostatic state at all relevant for the flow properties 
? It will be shown that while the stress fluctuations 
reflect the critical properties at point J, dynamical 
correlations are typical those observed in the flow 
of elasto-plastic solids. 

We will approach these questions starting with 
the quasistatic-flow regime. The advantage of qua- 
sistatic simulations is to provide a clean way of ac- 
cessing the transition region between elastic, solid- 
like behavior and the onset of flow. In the qua- 
sistatic regime flow is generated by a succession of 
(forcc-)equilibrated solid states. Thus, one can con- 
nect a liquid-like flow with the ensemble of solid 
states that are visited along the trajectory through 
phase-space. 

In Section l3~Tl we study the instantaneous statis- 
tical properties of the configurations generated by 
this flow trajectory at zero strain rate, and show 
that they display large fluctuations in several quan- 
tities, that are associated with the proximity to the 
jamming point. 

In Section 13.21 we follow the analysis of recent 
experiments^ and use a "four point correlation" 
tool to define a dynamical correlation length that 
characterizes the extension of the dynamical het- 
erogeneities observed in the flow process. This dy- 
namical length scale is shown to scale as the system 
size in the zero strain rate limit, independently of 
the distance to point J. The heterogeneity in the 
system is maximal for strains that correspond to 
the typical duration between the plastic avalanches 
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described above. 

We complement this analysis with preliminary 
results from dissipative molecular-dynamics sim- 
ulations that access strain-rates above the qua- 
sistatic regime. This allows us to assess the im- 
portance of dynamic effects in limiting access to 
certain regions of the landscape. Indeed, the re- 
sults at larger strain rate are system size indepen- 
dent, and reveal a surprising growth of the strength 
of the heterogeneities with increasing packing frac- 
tion away from <p c . 

2 Simulations 

Our system consists of N soft spherical particles 
with harmonic contact interactions 

E(r) = k(r-r c ) 2 . (1) 

Two particles, having radii fj and rj, only interact, 
when they are "in contact" , i.e. when their distance 
r is less than the interaction diameter r c = ri + rj . 
This system has been studied in several contexts, 
for example in. 16 i 26 i 32 The mixture consists of two 
types of particles (50 : 50) with radii r\ = 0.5d 
and r2 = 0.7d in two-dimensions. Three different 
system-sizes have been simulated with N = 900, 
1600 and 2500 particles, respectively. The unit 
of length is the diameter, d, of the smaller par- 
ticle, the unit of energy is kd 2 , where k is the 
spring constant of the interaction potential. We 
use quasistatic shear simulation, and compare some 
of the results with those obtained from dissipative 
molecular-dynamics simulations at zero tempera- 
ture. 

Quasistatic simulations consists of successively 
applying small steps of shear followed by a mini- 
mization of the total potential energy. The shear 
is implemented with Lee-Edwards boundary con- 
ditions with an elementary strain step of A7 = 
5 • I0~ 5 . After each change in boundary condi- 
tions the particles are moved affincly to define the 
starting configuration for the minimization, which 
is performed using conjugate gradient techniques.— 
The minimization is stopped when the nearest en- 
ergy minimum is found. Thus, as the energy land- 
scape evolves under shear the system always re- 
mains at a local energy minimum, characterized by 
a potential energy, a pressure p and a shear stress 
a. 



The molecular dynamics simulations were per- 
formed by integrating Newton's equations of mo- 
tion with elastic forces as deduced from Eq. ([1]) 
and dissipative forces 

Fij = -b [(vi - vj ) • fij] fij , (2) 

proportional to the velocity differences along the 
direction fij that connects the particle pair. The 
damping coefficient is chosen to be b = 1. Rough 
boundaries are used during the shear, the bound- 
aries being built by freezing some particles at the 
extreme ends in the y-direction, from a quenched 
liquid configuration at a given (f>. The system is 
sheared by driving one of the walls at a fixed ve- 
locity in the x direction, using periodic boundary 
conditions in this direction. 

For all system sizes, the distance between the 
top and bottom boundaries is 52.8c? and each of 
the boundaries has a thickness of A. 2d. The system 
size is changed by modifying the length of the box 
in the x-direction. 

3 Results 

3.1 Quasistatic simulations 

As is readily apparent from Fig.[TJ a typical feature 
of quasistatic stress-strain relations is the interplay 
of "elastic branches" and "plastic events". During 
elastic branches stress grows linearly with strain 
and the response is reversible. In plastic events the 
stress drops rapidly and energy is dissipated. 

In setting the elementary strain step, A7, care 
must be taken to properly resolve these events. Too 
large strain steps would make the simulations miss 
certain relaxation events. We chose a strain step 
small enough, such that most minimization steps 
do not involve any plastic relaxations. In conse- 
quence, the elastic branches are well resolved, each 
consisting of many individual strain steps. 

The succession of clastic and plastic events de- 
fines the flow of the material just above its yield- 
stress (j y ((/)) . The value of the yield-stress depends 
on volume-fraction and nominally vanishes at <f> c 
(see Fig. [5]). For finite systems, however, finite-size 
effects dominate close to 4> c such that one cannot 
observe a clear vanishing of a y . Rather, as Fig. [T] 
shows, one enters an intermittent regime, i.e. a fi- 
nite interval in volume-fraction in which the strcss- 
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Figure 1: Stress-strain relation for two different 
volume-fractions, <f> = 0.847 (top) and <j> = <\> c — 
0.8433 (bottom). At volume- fractions close to 4> c 
the signal is intermittent showing long quiet regions 
where the system flows without the building up of 
stress. 

signal shows a coexistence between jammed and 
"freely-flowing" states. This is evidence of a distri- 
bution of jamming thresholds, P{<fi c ), which sharp- 
ens with increasing the system-size. 16 i 28 A finite- 
size scaling analysis of this distribution allows one 
to extract the critical volume-fraction. To this end 
we count the number of jamming events that lead 
from the freely-flowing state to the jammed state 
and back 0. We find a maximum number of events 
at a certain (jj c (L), which can be extrapolated to 
L = oo to define the critical volume fraction of 
our simulation. The value we find, </> c = 0.8433 is 
slightly higher than what has been obtained pre- 
viously, however, evidence is mounting that 4> c is 
non-universal— and depends on the details of the 
ensemble preparation. Scaling properties in the 
vicinity of a jamming threshold, on the other hand, 
appear to be universal.— 

In Fig.[2]wc display the yield-stress a y , as a func- 
tion of volume- fraction, 5(f) = (f>—(f) c , and prcssurep. 
Finite-size effects are particularly strong when us- 
ing <j) as control variable. In the intermittent regime 
the average stress levels off to a system-size depen- 
dent value. 

For a closer illustration of a jamming event see the sup- 
porting material to our previous papcri^ 
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Figure 2: Average yield-stress as function of 
volume-fraction (left) and pressure (right). The 
yield-stress is determined as an average over stress- 
values just before plastic events occur, i.e at the 
top end of each elastic branch. 

Much better scaling behavior can be obtained, 
when using pressure as control variable, as this is 
characterized by the same finite-size effects as the 
shear-stress. In the following we will therefore use 
pressure as control variable. Be aware, however, 
that we do not run pressure-controlled simulations, 
as for example Peyneau and Roux— but use the 
average pressure, (p) ((/>), only to plot our simula- 
tion results. The value er/p « 0.1 obtained from 
Fig. [5] is consistent with these pressure-controlled 
simulations. On the other hand, the scaling with 
volume-fraction, p ~ S(j> , is slightly stronger than 
in linear elasticity at zero stress, 16 ' 20 where the 
pressure simply scales as 8<p. In view of the strong 
finite-size effects, the scaling with volume-fraction 
should, however, be taken with care. 

In the following we show results from five dif- 
ferent volume-fractions, 4> = 0.846, 0.848, 0.85, 
0.86 and <f> = 0.9, which are all above (j> c and 
outside the intermittent regime (i.e. no freely- 
flowing zero-stress states occur). For each volume- 
fraction we study three different system sizes with 
N = 900, 1600 and 2500 particles. 

3.1.1 Elastic properties 

As reviewed in the introduction a hallmark of the 
elasticity of solids in the vicinity of point J is the 
scaling of the linear elastic shear modulus g ~ 
p 1 / 2 . Similarly, the number of inter-particle con- 
tacts scale as z = zq + Ap 1 / 2 . 

We have analyzed the elastic branches in the 
steady-state flow to find (Figs. [3] and [4]) that the 
same scaling properties characterize the average 
nonlinear clastic modulus <7 aV g, which we define as 
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the local slope of the stress-strain curve, and also 
the associated contact numbers z avg . If we take 
these scaling properties as a signature of the criti- 
cality of point J, we can conclude that for the range 
of volume-fractions considered we are in the "criti- 
cal regime" . 
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Figure 3: (Top) Average nonlinear elastic shear 
modulus (?avg as function of pressure p. (Bottom) 
Probability distribution P{g) centered around av- 
erage value (? avg and rescaled width according to 
Ag = p°- 25 /N - 5 . Black solid line is a Gaussian 
pdf. 

As additional characterization of the ensemble of 
elastic states we report the probability distributions 
of shear moduli and contact numbers, respectively. 
Maybe surprisingly all the obtained distributions 
have approximately the same shape and can be su- 
perimposed on a single master curve. To achieve 
this we center each distribution around the average 
value and rescale the width with a factor p a N° . 

By looking carefully at the individual distribu- 
tions we do observe a slight trend towards the de- 
velopment of non-Gaussian tails close to <j) c . While 
non-Gaussian distributions are to be expected close 
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Figure 4: (Top) Average contact number z as func- 
tion of pressure p. The values for zq = z(p = 0) are 
determined from a best fit. They are smaller than 
z c = 4 as to the presence of rattlers, which have 
not been accounted for. (Bottom) Probability dis- 
tribution P(z) centered around average value z avg 
and rescaled width according to Az = p~ - 35 /TV ' 5 

to critical points,— the effect is quite small and all 
distributions have a well developed Gaussian core. 
The pronounced small-<? tail of P(g) is due to shear 
moduli that extend down to zero. Similar tails have 
been observed in— and related to a softening of the 
response upon approach towards plastic instabil- 
ities. Indeed, we found that manually suppress- 
ing states close to plastic events, the weight in the 
small-g tail is reduced. 

For the width of the g-distribution we obtain 
Ag = p0-25/7v - 5 . Thus, the absolute width of 
the distribution decreases with decreasing pressure, 
while the relative width, Ag/g avg diverges at point 
J. For the contact numbers, on the other hand, 
we find a divergence of the absolute width itself, 
Az = p~°- 35 /A^ - 5 . These enhanced fluctuations 
certainly support the view of Sz as an order pa- 



5 



rameter for a continuous jamming transition. The 
quantity Az 2 N would then be analogous to a sus- 
ceptibility, x ~ V ~ 7 > diverging with an exponent 
7 = 0.7. 

Our results arc different than those of Henkcs and 
Chakraborty,— where fluctuations of z are found to 
be independent of pressure, Az ~ p . Note, how- 
ever, the subtle difference in ensemble. These au- 
thors study a pressure-ensemble, in which states are 
generated by quenching random particle configura- 
tions to the local minimum of the potential energy 
landscape (similar to the procedure in Ref.— ). One 
may view these states as the inherent structures of a 
high temperature liquid. Our ensemble then corre- 
sponds to the inherent structures of a driven glassy 
material. We fix volume-fraction and sample only 
states that are connected by the trajectory of the 
system in phase-space. The ensemble therefore re- 
flects the dynamics of the system and the region of 
phase-space where it is guided to. 

3.1.2 Yield properties 

We now go beyond the properties of the elastic 
states and discuss aspects related to their failure 
during the plastic events. As indicated in the in- 
troduction, plastic events can be viewed as bifurca- 
tions in energy-landscape. A local energy minimum 
vanishes and the system has to search for a new 
minimum at lower energy and stress. In quasistatic 
dynamics this process is instantaneous. The asso- 
ciated stress-drop is therefore visible vertical 
line in the stress-strain relation (Fig. [T]). 

In the following we characterize the amount of 
dissipation during a plastic event by the associated 
stress-drop, Ac. The frequency of plastic events is 
discussed in terms of the length of elastic branches, 

A7avg ■ 

Just like the probability distributions within the 
clastic states, the functions P(Aa) and P(A7 avg ) 
(Figs. [5] and [5]) are universal and can be rescaled 
on a single master-curve. Here, it is sufficient to 
use the first moment of the distribution, i.e. the 
ensemble-averaged stress drops and elastic-branch 
lengths, respectively. As the logarithmic scale in 
the inset in Fig. [5] shows, the collapse for the stress- 
drop distribution is quite good for large as well 
as for small stress drops. The black line repre- 
sents a fit of the form Act -1 exp(— Act/cl) with 
the stress-scale ol ~ 5Aer avg . The intermediate 
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Figure 5: (Top) Distribution of stress-drops nor- 
malized with average values Aer avg . Inset shows the 
same figure in a log-log representation. (Bottom) 
Scaling of Aer avg with pressure p. 

power-law behaviour P ~ Act -1 reflects the lack of 
scale related to a typical event size. The only rele- 
vant scale is the exponential cut-off at <tl ■ Tewari 
et al— have reported an exponent of —0.7 in the 
energy-drop distribution at finite-strain rates. The 
simulated systems are somewhat smaller, however. 
Kabla et al.— have found an exponent of —1.5 in 
a vertex model for foams, in agreement with renor- 
malization group arguments.— 

The exponential tail has been observed in sev- 
eral different studies^ ' 11 ' 11 ' 41 ' 42 in two and in three 
spatial dimensions. Tsamados et al.— have further- 
more related this feature of the stress-drop distri- 
bution to the diversity of local flow-defects causing 
the plastic event. 

A similar universality has been observed by Mal- 
oney and Lemaitre.— Their simulations are con- 
ducted with three different interaction potentials 
but without changing the density, which is set to 
high values far away from the rigidity transition. 
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Figure 6: (Top) Distribution of elastic branch 
length normalized with average values A7 avg . (Bot- 
tom) Scaling of A7 avg with pressure p. The ratio 
g = Acr aV g/A7 avg , which defines a shear-modulus, 
is consistent with the scaling in Fig. [3] 



The authors have argued for a universal value of 
the "flow-strain" a y /g of a few percent. Appar- 
ently, this can only be true far away from <f> c . As 
the yield-stress vanishes faster than the shear mod- 
ulus, one finds a ratio a y /g ~ 5(f) 1 / 2 that vanishes 
at point J. Thus, particle configurations at the on- 
set of jamming are highly fragile and susceptible to 
even minute changes in the boundary conditions. 

The average stress-drop as well as the average 
length of elastic branches change with pressure and 
system-size as displayed in Figs. [5] and [6] @. As a 
function of pressure one observes an increase but 
with a slope that depends on system-size. The av- 



2 Note, that these values also depend on the elementary 
strain step used in the simulations. For any finite step-size 
there will be some small plastic events, that cannot be re- 
solved but only lead to an apparent reduction of the stress 
increase. This will artificially increase the length of clastic 
branches and decrease the weight of the small Atr-tail of the 
stress-drop distribution. 



erage stress-drops increase somewhat slower with 
pressure than the yield-stress H. The relative stress 
fluctuations Aa avg /a y are thus slightly enhanced at 
small pressures close to (j> c . The same trend is also 
visible in the total stress fluctuations as calculated 
by<(<7-(<7)) 2 >. 

The overall scale of both, Aer avg and A7 avg , de- 
creases with system-size to give a smooth stress- 
strain relation in the thermodynamic limit. Pre- 
vious studies ^ 10 i 11 have observed a scaling of the 
stress-drops with N~ 1,>2 . This includes^ a sys- 
tem of harmonically interacting particles, as stud- 
ied here, but at a rather high pressure. In gen- 
eral, we observe a weaker dependence on system- 
size, with an effective exponent that increases with 
pressure. Our data is consistent, however, with the 
value of 1/2 being the relevant high-pressure limit. 

3.2 Dynamical correlations 

Let us now turn to the dynamics of the system. In 
particular we want to characterize dynamic corre- 
lations in the motion of particles. While at volume- 
fractions above 4> c the isostaticity length-scale I* is 
clearly finite 0, there is nevertheless a large dynam- 
ical length-scale related to the flow arrest. This 
has, for example, been evidenced in a system of 
Lennard-Jones particles with dissipative dynam- 
ics.— We will show below that a similar length- 
scale occurs in our system of purely repulsively in- 
teracting particles, independent of the distance to 

4> c - 

Let us start by presenting the results from the 
quasistatic simulations. To define a dynamical cor- 
relation length we study heterogeneities in the par- 
ticle mobilities. To this end we use the overlap- 
function—^ 



/ 1 N 
(Q( 7 , a)) = ( — exp 

\ i=l 



"ma(7) 



2 1 



2a 2 



(3) 



of particles undergoing nonaffinc displacements 
Ui na during a strain interval of 7 0. Particles mov- 
ing farther than the distance a ("mobile"), have 



4 If we assume for the isostaticity length I 
would have values I* fss 5 at <b = 0.846 and I* P 



3 The effective exponents range from 0.8 to 0.9 as com- 
pared to an exponent 0.95 for the yield-stress 

= 1/Sz, wc 
1 1 at (j> = 0.9 
(^-values taken from Fig. [4}. 

5 To calculate the overlap function we use the non-affinc 
displacements in the gradient direction (irrespective of the 
distance to the wall). 
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Q m 0, while those that stay within this distance 
("immobile") have Q ~ 1. 

As a function of strain 7, the average over- 
lap (Q) will decay, when particle displacements 
u na are comparable to the probing length-scale a. 
The overlap function is similar to the intermedi- 
ate scattering function with wave- vector q ~ 1/a. 
Thus, a sets the probing length-scale. The decay 
of Q(j,a) then gives an associated structural re- 
laxation strain, 7* (a), on which particle positions 
decorrelate. 

In the following we are interested in the dynami- 
cal heterogeneity of Q and the fluctuations around 
its average value 

X 4(a, j) = N (<Q( 7 , a) 2 ) - (Q a ( 7 , a)) 2 ) , (4) 

which defines the (self-part of the) dynamical sus- 
ceptibility Xi- This is displayed in Fig. [7] func- 
tion of both strain 7 and probing length-scale a. 
For each 7 it has a well defined peak (at a* (7)) that 
parallels the decay of the overlap function (Q)-— 
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Figure 7: Dynamical susceptibility \4 as function of 
probing length-scale a for various different strains 
7 = (5,20,50,200,500,2000) • 10~ 5 (from left to 
right) and <j> = 0.9. The maxima of the curves 
(black circles) define the amplitude h(j). 

The strength of the correlations are encoded in 
the peak-height, ^,(7) = Xi{ a * il) 1 1) (black circles 
in Fig. [7])). As Xi can be written as the integral 
over a correlation function, it is connected to the 
correlation volume, or to the number of correlated 
particles. Assuming that this volume forms a com- 
pact region in spac e 45 ' 46 we can relate the am- 
plitude of Xi to a dynamic correlation length via 
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Figure 8: Amplitude h(j) as taken from the qua- 
sistatic simulations. Data for different volume- 
fractions and system-sizes. The axes are nor- 
malized according to the scaling form h(p/) = 
A r /i(7/A7 avg ), with A7 avg taken from FigHJl 

Following the maxima in Fig.[7Jfrom left to right, 
one sees that the amplitude first increases and then 
quickly drops to small values. This implies that 
there is a finite strain 7, at which Xi presents an ab- 
solute maximum. To extract this maximum we plot 
in Fig. [8] the amplitude ^1(7) for various volume- 
fractions (f> and system-sizes N. 

There are two surprising features in this plot. 

First, by rescaling the strain-axis with the aver- 
age length of elastic branches, A7 avg (see Fig. [5]) 
we find reasonable scaling collapse for all studied 
volume-fractions and system-sizes. This implies 
that cooperativity, as measured by the amplitude of 
Xi and the length of clastic branches are intimately 
related. The frequency of plastic events sets the 
strain-scale for dynamical heterogeneities. As the 
length of the elastic branches decreases with system 
size, the absolute maximum shifts towards smaller 
strains, with ^1(7) becoming effectively a decreasing 
function of 7 in the thermodynamic limit. 

The second surprising feature in Fig. [5] is the 
system-size dependence of h. It turns out that 
h/N rather than h itself is independent of system- 
size, indicating a finite variance of the distribu- 
tion of Q values in the thermodynamic limit (see 
Eq. Q). Assuming the connection with the corre- 
lation length to hold, £ 2 ~ h, this implies a cor- 
relation length that is proportional to the length 
of the simulation box, £ w 0.3L, independent of 
volume-fraction and distance to point J. This illus- 
trates the fact that quasistatic dynamics is inhcr- 
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Figure 9: Xi(j) f° r different strain-rates 7 and 
a = 0.01. comparison with quasistatic ('qs') sim- 
ulations. Note the different boundary conditions 
used. MD simulations arc with walls, while qua- 
sistatic simulations have periodic boundary condi- 
tions. This may explain the difference in the am- 
plitude. 



ently dominated by system-size effects, as already 
shown in previous works . 14 i 15 i 28 i 47 Note, that this 
system-size dependence is of different origin than 
the finite-size effects present within the above men- 
tioned intermittent regime, which occurs close to 
4> c . The intermittency can be avoided by staying 
away from <f> c . In contrast, the system-size de- 
pendence encountered here is quite independent of 
volume-fraction, but rather a generic feature of the 
quasistatic regime, as we show now. 

To this end let us turn to the molecular-dynamics 
simulations. We will show that the dependence 
on system-size indeed reflects the saturation of a 
length-scale that is finite for larger strain-rates and 
increases towards the quasistatic regime . As 
Fig. [9] shows, the amplitude of \i increases when 
reducing the strain-rate (a = 0.01, <f> = 0.9) and ap- 
proaches the quasistatic limit for small strain-rates. 
Also the strain 7 m , at which Xi is maximal is very 
well reproduced in the dynamic simulation. 

Fig. [10] demonstrates a saturation of the ampli- 
tude Xm = Xiilm) at small strain-rates indicating 
that the quasistatic regime is entered. Comparing 
with the quasistatic simulation, we find somewhat 
smaller values for the amplitude. It should be re- 
membered, however, that different boundary con- 



6 A more detailed account of these simulations will be 
presented in: P. Chaudhuri and L. Bocquct, in preparation 
(2010). 
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Figure 10: (Left) Peak-height Xm = X^ilm) as 
determined from Fig. [9] The saturation at small 
strain-rates is an indication of the quasistatic limit , 
in which Xm ~ N. In contrast, at high strain- 
rates no significant dependence on system-size is 
observed. (Right) Peak-height Xm as function of 
volume- fraction. 



ditions have been used. The rough walls used in 
the molecular dynamics simulations arc likely re- 
sponsible for the reduction of the peak-height as 
compared to the quasistatic simulations (which are 
performed with periodic boundary conditions). 

The presence of the quasistatic regime is also ev- 
idenced by the fact that the amplitude Xm within 
the plateau depends on system-size, just as in the 
quasistatic simulations. Outside this regime, on 
the other hand, no significant A-dependence is ob- 
served. In effect this means that the quasistatic 
regime shrinks with increasing system-size. The 
strainrate 7 qs (A) that describes the crossover to 
the quasistatic regime decreases with N . This is 
in line with Refs., 14 i 15 where a power- law depen- 
dence 7 qs ~ 1/A is reported. From our data we 
cannot make any definitive statement about this 
dependence. 

The results arc furthermore consistent with those 
of Ono et ai— The lowest strain-rate accessible in 
the latter study was 7 = 0.0001. At this strain- 
rate the correlation length was observed to be on 
the order of 3 in agreement with our data. 

We also probed the volume-fraction dependence, 
by performing runs at (j> = 0.85,0.87,0.9,0.95 and 
1. The resulting amplitude of Xi is given in Fig. llQI 
Interestingly we observe a mild increase in the am- 
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plitude with volume-fraction, signalling enhanced 
correlations away from <f> c . 

This is not due to the special choice of the pa- 
rameter a = 0.01. We have found the same trend 
when fixing 7 and viewing Xi as function of a as 
for the quasistatic simulations in Fig. [7] Finally, we 
have also calculated the absolute maximum of xa, 
viewed as function of both a and 7. In all cases, the 
amplitude increases with volume-fraction. Both, 
the increase of the length-scale with lowering the 
strain-rate and the increase with volume-fraction 
are consistent with the recently proposed elasto- 
plastic model of Bocquet et al— 

Given the trend in Fig. 1101 one may specu- 
late about a vanishing dynamical correlation length 
(taken at constant strain-rate 7 = 10~ 4 outside the 
quasistatic regime), as </> c is approached. Such a 
behavior is indeed compatible with our data and 
has recently been observed in the rheology of a 
concentrated emulsion confined in gaps of differ- 
ent thickness.— Here we interpret this surprising 
feature in the following way: at a given packing 
fraction, dynamical correlations increase with de- 
creasing 7, and saturate at a TV-dependent value in 
the quasistatic regime. This cross-over to the qua- 
sistatic regime does not only depend on system-size 
but also on packing fraction: 7 qs (TV, </>) . 

For (j) closer to <j) c the energy landscape becomes 
increasingly flat. Particle relaxations take longer— 
and smaller strain-rates arc needed to allow for full 
relaxation into the local energy minimum. Thus, 
smaller strain-rates are needed to reach the qua- 
sistatic behavior and 7 qs decreases towards <p c . 
Hence, reducing <f> towards 4> c at a fixed strain rate 
is " equivalent" to increasing the strain rate relative 
to 7 qs , and results in a decrease of dynamical corre- 
lations. In contrast, corrrelations taken at a strain 
rate 7 = 7 qs (0), are independent of <p and given by 
their quasi-static values, as displayed in Fig. [5] 

4 Discussion and Conclusion 

We have discussed the small strain-rate elasto- 
plastic flow of an athermal model system of soft 
harmonic spheres. In particular, we were inter- 
ested in the flow properties at and above a critical 
volume- fraction (point J), at which the yield-stress 
of the material vanishes. This regime combines the 
more traditional elasto-plastic flow of solids above 



their yield-stress with the breakdown of the rigidity 
of the solid state at point J. 

We found that this breakdown is visible in the en- 
semble of states visited during a flow simulation in 
a similar way as in the linear elasticity of the solid. 
For example (Fig. @|, we showed that the average 
number of particle contacts scale with the square- 
root of pressure, just as in linear elasticity In con- 
trast, the fluctuations around this average value 
show a distinct behavior that has not been observed 
previously. We showed that the contact-number 
fluctuations actually diverge upon approaching the 
critical volume-fraction from above, making the 
contact number an ideal candidate for an order 
parameter of a continuous jamming transition as 
observed under steady shear. The relative fluctua- 
tions of the shear modulus and those of the shear 
stress also diverge in the same limit. Going beyond 
the characterization of the average elastic proper- 
ties we have studied the statistics of plastic events 
(Figs. [5] and [6]). It seems that all distributions have 
universal scaling forms reminiscent of standard crit- 
ical phenomena. 

From all these results, it would be tempting to 
say that it is the energy landscape whole 
that becomes critical at point J. Isostatic elastic- 
ity would then be just one aspect of this criticality, 
another one could be the intermediate power-law 
tail in the stress-drop distribution. This critical as- 
pect is also illustrated by the intcrmittency in the 
stress response of finite-size systems (see Fig|T]) and 
by the growth of an isostatic correlation length in 
the quasistatic response when point J is approached 
from below.— 

At strain-rates above the quasistatic regime, the 
dynamics limits access to certain regions of the en- 
ergy landscape. While the dynamics is still highly 
correlated, the dynamical correlation length, as 
measured by the amplitude of the four-point sus- 
ceptibility X4y remains finite and actually decreases 
with lowering the volume- fraction towards 4> c . 

In the quasistatic regime we have shown that Xi 
reflects, in two ways, the interplay of elastic load- 
ing and plastic energy release (Fig. [5]). First, the 
typical strain-scale of heterogeneity is set by the fre- 
quency of plastic events. Second, the amplitude of 
X4 scales with system-size, which highlights the fact 
that the quasistatic, plastic flow regime is, in fact, 
a finite-size dominated regime with a correlation 
length that is limited by system size. This behavior 
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should be contrasted with the one observed below 
4> c , where a large but finite correlation length has 
been identified, which is governed by the approach 
to point J.— 

Upon increasing the strain-rate we have shown 
that the correlation length starts to decrease out- 
side the finite-size scaling regime (Fig. [TU|) . Olsson 
and Teitel— infer from their flow simulations that 
shear-stress should be viewed as a "relevant pertur- 
bation" to point J, such that a different fixed-point 
and indeed different physics is relevant for the flow 
behaviour at finite stress. Our findings support this 
picture for the dynamical correlations, which ap- 
pear to behave similarly to those observed in mod- 
els of elasto-plastic flo w 13 i 14 i 52 or in low temper- 
ature glasses: 10 i 15 correlations increase upon low- 
ering the strain-rate and saturate at a system-size 
dependent value in the quasistatic regime. 

The flow behaviour in the vicinity of point J is 
therefore influenced by a complex combination of 
two critical behaviour. Large stress fluctuations 
(relative to the yield stress) , or geometrical changes 
(number of neighbours) reflect the enhanced sen- 
sitivity of the material to small changes in exter- 
nal conditions at point J, and are specific proper- 
ties of the energy landscape at this point. On the 
other hand, dynamical correlations above point J 
are dominated by the system size, and build up 
progressively as the strain rate is decreased, as in 
any elasto-plastic system, and are not particularly 
sensitive to the proximity of point J. 
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